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ABSTRACT 


The  Pseudo  Wigner- Ville  Distribution  is  a  signal  transformation  of  an  input  time  signal 
into  a  joint  time-frequency  domain,  that  provides  an  excellent  characterization  of  an  input 
signal  as  well  as  it's  respective  energy  content.  A  smoothing  window  and  energy  sensitivity 
analysis  of  the  discretized  Pseudo  Wigner- Ville  Distribution  is  presented  along  with  a 
machinery  monitoring  application  using  the  resulting  Wigner- Ville  energy  of  the  sampled 
signal.  The  ability  to  apply  the  Pseudo  Wigner-Ville  Distribution  to  both  steady-state  and 
transient  machinery  may  enable  the  monitoring  and  diagnostics  of  virtually  any  piece  of 
equipment. 
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I.     INTRODUCTION 

The  recent  world  events  have  caused  drastic  changes  in  many  areas  and  made  the  future 
that  much  more  unpredictable.  With  this  high  degree  of  uncertainty  and  drastic  financial 
cutbacks  as  evidenced  in  the  Defense  budget,  the  emphasis  of  technological  advance  and 
dominance  is  all  the  more  heightened.  With  the  requirement  to  "do  more  with  less",  it  is 
imperative  to  find  ways  to  save  money.  It  is  with  this  mind-set  that  better  ways  to  prolong 
the  life  of  operating  machinery  and  defray  unneeded  maintenance  costs  be  investigated.  For 
many  years,  most  transient  machinery  has  gone  unmonitored  and  removed  at  periodic  time 
intervals  for  refurbishment  regardless  of  the  current  machinery  condition.  Developing  a 
method  to  monitor  and  diagnose  these  machines  of  a  transient  nature  will  provide  a  way  to 
determine  when  overhaul  or  replacement  is  needed. 

Vibrational  data  has  long  been  used  as  the  tool  to  measure  a  machine's  health.  Early 
methods  of  monitoring  machinery  had  been  for  the  experienced  personnel  to  listen  or  touch 
the  machine  to  detect  if  a  fault  exists,  however,  this  method  is  dependent  on  human  senses 
and  the  same  individual  conducting  the  monitoring.  It  was  for  this  reason  that  vibration  levels 
as  monitored  by  a  dynamic  signal-analyzer  are  now  used  for  machinery  health  monitoring. 
The  vibration  levels  of  a  good  condition  establish  a  baseline  to  allow  a  means  of  comparison 
when  subsequent  measurements  are  taken.  These  dynamic  signal  analyzers  are  time 
independent  and  do  not  show  the  time  dependency  required  by  transient  machinery.  It  is 
because  of  this  limitation,  that  a  three-dimensional  time,  frequency  and  energy  representation 
such  as  the  Pseudo  Wigner-Ville  Distribution  is  being  investigated  as  a  method  to  monitor 
transient  machinery. 


A.  PREVIOUS  WORK  AND  ACCOMPLISHMENTS 

This  thesis  is  a  follow-on  work  which  has  been  preceded  by  LT  G.  Rossano  and  LCDR 
J.  Calahorrano,  both  in  conjunction  with  Prof.  Y.S.  Shin.  Rossano  established  the  initial  form 
of  the  computer  program  that  calculates  the  discretized  Pseudo  Wigner-Ville  Distribution. 
He  evaluated  the  effect  that  different  signals  and  procedures  had  on  the  Distribution  along 
with  characterizing  an  actual  pump  signal.  Calahorrano  made  some  further  improvements  to 
the  computer  program  and  then  applied  this  improved  version  of  the  computer  program  to 
characterize  some  classified  signals.  Both  of  these  works  have  been  a  great  benefit  in  dealing 
with  characterizing  signals,  but  have  not  investigated  the  smoothing  method  and  energy  or 
amplitude  of  the  discretized  Pseudo  Wigner-Ville  Distribution. 

B.  OBJECTIVES  OF  CURRENT  RESEARCH 

The  Pseudo  Wigner-Ville  Distribution  provides  an  outstanding  time-frequency  domain 
spectra  representation  of  an  input  signal  as  shown  through  the  works  among  others  of  Rossano 
and  Calahorrano  here  at  the  Naval  Postgraduate  School.  The  purpose  of  this  research  is: 

•  To  review  and  modify  the  computer  program  developed  here  at  the  Naval  Postgraduate 
School  to  support  the  smoothing  window  and  energy  analysis  work. 

•  To  conduct  a  smoothing  window  sensitivity  analysis. 

•  To  investigate  the  energy  content  of  the  discretized  Pseudo  Wigner-Ville  Distribution. 

•  To  apply  the  developed  smoothing  window  and  energy  principles  to  an  actual 
machinery  monitoring  application. 

•  To  investigate  the  feasibility  of  using  the  Pseudo  Wigner-Ville  Distribution  energy  as 
a  pre-processed  input  to  a  Neural  Network. 


II.     THE  PSEUDO  WIGNER-VILLE  DISTRIBUTION 

A.       THEORETICAL  BACKGROUND 
1.       Evolution 

The  current  day  Pseudo  Wigner-Ville  Distribution  is  a  three  dimensional  (time, 
frequency,  and  energy)  representation  of  a  signal  that  is  particularly  well  suited  for  analysis 
of  non-stationary  signals.  Eugene  Wigner  [Ref.  1]  first  introduced  the  Wigner  distribution  in 
1932  to  study  the  problem  of  statistical  equilibrium  in  the  area  of  quantum  mechanics.  His 
work  was  furthered  in  1948  by  J.  Ville  [Ref.  2],  who  used  the  Wigner  distribution  in  the  area 
of  signal  analysis.  A  major  contribution  of  Ville's,  is  the  use  of  analytic  signals  as  an  input 
to  the  distribution  vice  the  customary  real  signal.  The  advantages  of  using  an  analytic  signal 
in  the  distribution  are  two-fold  [Ref.  3].  First,  the  distribution  of  a  real  signal  results  in  a 
symmetrical  spectrum  with  only  half  of  the  result  containing  useful  information,  while  the 
use  of  an  analytic  signal  avoids  this  negative  frequency  redundancy.  Secondly,  by  accounting 
for  only  the  positive  frequencies  it  satisfies  both  the  practical  and  the  mathematical 
completeness  of  the  problem.  The  third  part  of  the  Pseudo  Wigner-Ville  Distribution  (ie. 
Pseudo)  ,  arises  from  the  need  to  apply  a  smoothing  window  to  the  resulting  distribution. 
This  is  a  result  of  the  presence  of  cross  terms  that  arise  from  multiple  component  signals  and 
will  be  covered  in  depth  later. 

The  works  of  T.A.C.M.  Claasen  and  W.F.G.  Mecklenbrauker  [Refs.  4,5,  and  6]  in 
1980  have  paved  the  way  for  many  recent  applications  of  the  Pseudo  Wigner-Ville 
Distribution.  Their  three  part  paper  is  an  all  encompassing  work  that  has  served  to  highlight 
the  capabilities  of  the  Pseudo  Wigner-Ville  Distribution  and  allow  for  greater  knowledge  of 
the  subject  with  an  emphasis  on  signal  processing  and  analysis. 


A  considerable  amount  of  recent  work  using  the  Pseudo  Wigner- Ville  Distribution 
have  come  in  the  areas  of  optics  [Refs.  7,8,  and  9]  and  speech  [Ref.  10  and  1 1].  Additionally, 
T.J.  Wahl  and  J.S.  Bolton  of  Purdue  University  [Ref.  12]  have  used  the  Pseudo  Wigner- Ville 
Distribution  to  analyze  structural  impulse  responses.  Two  recent  papers  have  investigated  the 
use  of  the  Wigner- Ville  Distribution  in  the  area  of  machinery  monitoring.  Flandrin  et.  al. 
[Ref.  13]  proposed  the  Wigner- Ville  as  a  means  to  confirm  a  machinery  diagnosis  in  the 
specific  application  of  a  Pressurized  Water  Reactor  power  plant  and  lastly,  Forrester  [Ref.  14] 
has  investigated  the  Wigner- Ville  Distribution  as  a  method  for  fault  detection  in  gears.  As 
evidenced  by  the  cited  examples,  the  capability  and  adaptability  of  the  Wigner-Ville 
Distribution  is  enormous  with  applications  of  it's  use  rapidly  expanding. 
2.       Function  Definition 

The  Wigner-Ville  Distribution  is  a  transformation  of  a  signal  into  the  time- 
frequency  domain.   By  definition,  the  cross  Wigner-Ville  Distribution  is  defined  as: 


WDFzs-fe-i™  ^(t+-|)    s*(t--|)    dx 


where:         r(t)  =  a  complex  time  history 

s(t)  =  a  complex  time  history 

t  =  time 

w  =  frequency 

*  =  complex  conjugate 
Similarly,  the  auto  Wigner-Ville  Distribution  is  defined  as: 

oo 

WDFr-    /VJut  r(t+  — )   r*(t-—)   dx 
r'r       J  2  2 


Since  the  concentration  of  this  research  is  with  the  representation  of  a  single  input  signal,  the 
auto  Wigner- Ville  Distribution  will  be  used.  This  expression  is  essentially  a  Fourier  transform 
of  the  auto-correlation  of  a  signal,  which  may  be  thought  of  as  a  three  dimensional  spectral 
density  function. 

3.       Distribution  Properties 

There  are  several  properties  of  the  Distribution  that  are  worth  noting.  The 
properties  listed  below  have  been  shown  in  many  works  concerning  the  Wigner- Ville 
Distribution.  These  properties  will  be  shown  in  the  form  of  the  auto  Wigner-Ville 
Distribution  [Ref.  15]. 

Time  and  frequency  shift  properties  follow: 

•  A  time  shift  in  the  signal  corresponds  to  a  time  shift  of  the  Distribution: 

WD^(t-T)S(t-x)(t/CO)    -  mFe§B{t-x,v) 

•  A  frequency  shift  in  the  signal  corresponds  to  a  frequency  shift  in  the  Distribution: 

WDF6jq:      j0t3(ttu)    -   WDF      (t,Q-Q) 


•    Both  a  frequency  and  time  shift  in  the  signal  corresponds  to  the  same  in  the 
Distribution: 


W^e^sU-tKe^sU-t)^'")     "    WDFs.  s  <  t"T  ,  G>-G  \ 


Since  our  concern  is  in  the  area  of  machinery  monitoring  and  diagnostics,  the  aforementioned 
properties  are  a  nessesity  in  the  development  of  a  suitable  monitoring  program. 


The  next  three  properties  provide  information  concerning  the  energy  contained 

within  the  Wigner-Ville  Distribution. 

•    The    integration    of   the    Distribution    over   the    frequency    domain    provides    the 
instantaneous  signal  power: 


-^-  fwDFs  s(t(w)dw  -ls(t)   P 
2tc  J  s,s 


The  integration  of  the  Distribution  over  the  time  domain  provides  the  energy  density 
spectrum: 


fwDF8tS(t,w)  dt  -  I  S(«)    I2 


•    The  integration  over  both  the  time  and  frequency  domain  provides  the  total  energy  of 
the  input  signal: 


—  [  (wdf^  At,u>)  dtdu  -Ms 
2ti  J  J        s,s 
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Again,  as  with  the  time  and  frequency  properties,  the  energy  contained  in  the  Wigner-Ville 
Distribution  must  be  definable  in  order  to  use  this  in  a  machinery  monitoring  program. 
Additional  properties  may  be  found  in  Reference  4. 
4.       Cross  Terms  Within  The  Distribution 

A  particularly  troublesome  characteristic  of  the  Wigner-Ville  Distribution  is  the 
presence  of  cross  terms.  These  cross  terms  result  in  the  Distribution  taking  on  negative  values 
as  shown  later  in  Figure  2.  Since  the  properties  above  have  shown  that  the  Wigner-Ville 
Distribution  represents  an  energy  value,  these  negative  values  raise  an  area  of  concern.  Jones 


and  Parks  [Ref.  16]  show  that  these  cross  terms  are  essentially  the  interference  of  multiple 
component  signals  that  occur  when  there  is  a  non-zero  cross-correlation  within  the  signal. 
The  equation  listed  below  shows  the  cross-term  that  results  from  a  signal  that  consists  of  two 
components. 

WDFx+y(  t,0>)    -   WDFx(t,U>)    +  2RE[WDFXty(tfa))    +   WDFy{t,u) 

As  shown,  the  cross-term  is  represented  by  the  center  term  while  the  other  terns  represent 
their  respective  signal  components.  There  are  several  methods  for  eliminating  these  terms  and 
providing  a  more  presentable  Distribution  and  there  has  been  a  considerable  amount  of 
research  in  the  are  of  an  optimum  smoothing  method  for  the  Wigner-Ville  Distribution  as 
indicated  by  References  3,4,17,  and  18.  The  method  chosen  for  smoothing  the  Wigner-Ville 
Distribution  in  this  work  is  a  post  processing,  sliding  Gaussian  window  function.  What  is 
meant  by  post  processing,  is  that  the  distribution  is  smoothed  after  calculation  of  the  Wigner- 
Ville  Distribution. 

5.       Gaussian  Smoothing  Function 

The  Gaussian  window  function  was  first  proposed  mathematically  as  a  means  to 
smooth  the  Wigner-Ville  Distribution  by  N.D.  Cartwright  [Ref.  19].  The  continuous  form  of 
the  Gaussian  window  function  follows: 

G(t,u) i e     2o*     2al 

Wahl  and  Bolton  in  Reference  12  have  furtherd  this  work  and  developed  the  sampled  form 
of  the  window  function  that  is  used  in  our  discretized  form  of  the  Pseudo  Wigner-Ville 


Distribution  for  smoothing.  The  effect  of  various  size  smoothing  windows  will  be  presented 
in  greater  depth  in  Chapter  III. 

B.        WIGNER-VILLE  DISTRIBUTION  COMPUTER  PROGRAM 

The  computer  program  listed  in  the  Appendix  is  a  method  that  was  developed  here  at 
the  Naval  Postgraduate  School  for  calculating  the  discretized  Pseudo  Wigner-Ville 
Distribution.  As  mentioned  earlier,  the  program  that  computes  the  discretized  Pseudo 
Wigner-Ville  Distribution  is  largely  the  work  of  G.  Rossano  [Ref.  15],  who  initially 
investigated  the  use  of  the  Wigner-Ville  Distribution  as  a  tool  for  machinery  monitoring.  For 
the  purposes  of  this  work,  the  program  has  been  thoroughly  reviewed  and  rewritten  to  allow 
for  ease  of  readibility  and  modification.  The  current  version  enables  the  user  to  look  at  not 
only  the  resulting  smoothed  version  of  the  Pseudo  Distribution,  but  also  the  unsmoothed 
Wigner-Ville  Distribution.  The  second  major  modification  was  to  program  the  smoothing 
window  subroutine  to  allow  for  the  use  of  different  size  windows  in  smoothing  the 
Distribution.  These  changes  have  proven  invaluable  for  evaluating  the  effect  of  various  sized 
smoothing  windows  and  in  conducting  the  energy  sensitivity  analysis.  Presented  below  are 
the  two  major  equations  that  comprise  the  calculation  of  the  discretized  Pseudo  Wigner-Ville 
Distribution  as  performed  in  the  program  listed  in  the  Appendix. 

1.       Discrete  Time  Wigner-Ville  Distribution 

As  mentioned  earlier,  the  Wigner-Ville  Distribution  is  essentially  a  Fourier 
transform  of  an  auto-correlation  function.  There  have  been  a  couple  of  equations  developed 
to  represent  the  discrete  Wigner-Ville  Distribution.  The  version  that  has  been  used  in  the 
development  of  the  program  listed  in  the  Appendix,  is  one  that  was  developed  in  Reference 
4  by  Classen  and  Mecklenbrauker  and  and  later  used  in  the  work  of  Wahl  and  Bolton  [Ref. 
12]  and  Rossano  [Ref.  15].  The  only  difference  of  note  between  the  various  discrete  versions 
is  essentially  a  difference  in  constants.  All  of  the  versions  are  basically  a  Fourier  transform 


of  an  auto-correlation  function.  Shown  below  is  the  discretized  equation  that  has  been  used 
for  programming  purposes  in  this  work: 

N 

WDFss{itnx>,jM)    -  J^i?e(FFr[2At  CORR(i)]) 

j'-i 

where:         Re  =  Real  part 

Corr(i)  =  complex  auto-correlation  of  the  signal  s(t) 
FFT  =  Fast  Fourier  Transform 
At  =  sample  time 

2.       Discretized  Gaussian  Window  Function 

Due  to  the  need  to  deemphasis  the  resulting  cross-terms  within  the  Wigner-Ville 
Distribution,  a  suitable  method  for  smoothing  has  been  obtained.  The  sampled  equation  of 
the  selected  smoothing  window  is  shown  below  [Ref.  12]: 


G(t,  o>)    -  = e 

2rcoto-w 


where:         -2M<m<2M       &       -2N<n<2N 


ot  =  MAt       &       ou  =  NAw 


(mAt)2    +     (nAu)2  , 
2o2r  2ai 


III.     SMOOTHING  WINDOW  ANALYSIS 

As  shown  above,  the  smoothing  function  or  window  that  is  used  in  this  work  is  a 
Gaussian  window  that  is  applied  to  the  unsmoothed  Wigner-  Ville  Distribution.  Figure  1  is  the 
raw  time  input  signal  of  a  100  &  400  Hz  sine  wave  that  will  be  processed  through  the  Wigner- 
Ville  Distribution.  Figure  2  is  the  unsmoothed  Distribution  that  results  from  this  combined 
100  and  400  Hz  sine  wave.  As  can  easily  be  seen  in  Figure  2,  the  resulting  cross  terms  are 
very  dominant  and  need  to  be  deemphasized.  This  Gaussian  smoothing  function  is  applied 
after  the  Wigner-Ville  Distribution  has  been  calculated  (ie.  post  processing),  which  allows  for 
visual  inspection  of  both  the  unsmoothed  Wigner-Ville  Distribution  and  the  effect  that 
different  sized  windows  may  produce  on  the  smoothed  Distribution  (ie.  Pseudo).  For 
example,  Figure  3  is  the  smoothed  result  of  Figure  2  after  applying  a  13  by  13  sized  window. 
The  next  few  sections  will  discuss  some  of  the  characteristics  of  this  post-processing 
smoothing  window  and  what  may  be  gained  or  lost  through  varying  the  size  of  the  applied 
window. 
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Figure  1.    100  &  400  Hz  Sine  Wave  Input 
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100  *  400  Hz  SINE  WAVE 
UNSMOOTHEP  DISTRIBUTION 
AMPLITUDES  =  1.0~ 


Figure  2.   Unsmoothed  Wigner-Ville  Distribution 
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1Q0  *  400  Hz  SINE  WAVE 

SMOOTHING  WINDOW  ;  la^ig 

AMPLITUDES  =  1.0 


Figure  3.     Pseudo  Wigner-Ville  Distribution  (13  X  13  Window) 
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A.       SIZE  CRITERION 

As  shown  by  N.D.  Cartwright  [Ref.  19]  for  the  continuous  form  of  the  Wigner-Ville 
Distribution,  the  size  of  the  smoothing  window  may  be  selected  such  that  a  positive 
Distribution  will  result.  Following,  are  these  criterion  that  have  been  mathematically  proven 
for  the  continuous  case  to  assure  a  positive  distribution: 

where: 

o,.-  MA  t  o„  -  itfAa> 

By  chosing  M  and  N  such  that  this  criterion  is  met,  the  resulting  smoothed  discretized 
Wigner-Ville  Distribution  does  not  necessarily  result  in  all  of  the  values  being  positive.  This 
may  be  seen  in  Figure  4,  where  the  window  that  was  used  was  again  a  13  by  13  (ie.  N  =  13 
&  M  =  13)  sized  window.  However,  it  does  deemphasize  the  large  cross-terms  to  the  point 
that  they  are  nearly  zero  and  yields  a  very  presentable  Pseudo  Wigner-Ville  Distribution.  The 
ability  to  deemphasize  these  cross  terms  is  a  necessity  if  the  desire  is  to  evalute  the  energy 
change  within  a  signal.  For  the  remainder  of  this  work,  the  window  size  selected  (ie.  N  &  M, 
N  X  M,  or  N  by  M)  will  be  such  that  this  criterion  is  met  as  closely  as  possible. 
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100  &  400  Hz  SINE  WAVE 

SMOOTHING  WINDOW  ;  13  X  13 

AMPLITUDES  =  1.0 


Figure  4.    Smoothed  Distribution  with  Negative  Values 
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B.        TIME  AND  FREQUENCY  DEPENDANCE 

Since  the  criterion  shown  previously  are  time  and  frequency  dependent,  there  must  be 
a  relation  between  the  two  that  can  be  reduced  to  one  parameter.  The  customary  time  and 
frequency  relations  of  signal  analysis  are  slightly  changed  in  the  context  of  the  Wigner- Ville 
Distribution.  Shown  below  are  the  time  and  frequency  relations  when  working  with  the 
Wigner- Ville  Distribution. 

-  DPxAt  ;  Af - ;  fM„  -  2xDPxAf 


-max  -  4xt 


max 


where:         tmax  =  maximum  signal  length       &       fmax  =  maximum  frequency 
At  =  sampling  time       &       Af  =  frequency  resolution 
DP  =  number  of  data  points 

The  definitions  of  the  size  criterion  discussed  earlier  for  the  smoothing  window  follow: 


2    2 


and  since: 


a\  -    (MM)2  oi  -    (JVAg>)2 

and  the  time  and  frequency  are  related  in  the  Wigner-Ville  Distribution  through  the  number 
of  data  points  as  shown: 

Ao) 5L-j- 

2  DP  At 
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Making  this  substitution  into  the  size  criterion,  yields  a  very  simple  result  for  determining  the 
size  of  the  smoothing  window  required. 

MxN  *   — 

In  looking  at  this  result  it  seems  as  if  the  window  is  completely  independent  of  time  and 
frequency,  but  it  must  be  rememberd  that  M  and  N  are  parameters  that  specify  the  lengths 
of  the  smoothing  window  in  the  time  and  frequency  directions  respectively. 

C.        RELATION  TO  NUMBER  OF  DATA  POINTS 

Another  area  of  interest  is  the  relation  of  the  window  size  to  the  number  of  data  points 
in  the  sampled  signal.  The  best  way  to  show  this  relation  is  through  an  example.  For 
instance,  a  signal  of  a  length  4.0  seconds  that  is  sampled  using  512  data  points,  will  result  in 
a  At  =  7.8125  msec  and  Af  =  0.0625  Hz  as  calculated  using  the  constraints  of  time  and 
frequency  as  defined  by  the  Wigner-Ville  Distribution.  Since  the  example  has  used  512  data 
points,  the  previous  section  has  shown  that  chosing  N  =  10  and  M  =  16  will  approximately 
satisfy  the  size  criterion  limitations.  Additionaly,  by  definition  of  the  smoothing  window,  the 
area  of  coverage  will  be: 

-2M±m<L2M         -2N<.n<L2N 

This  means  that  for  the  example  shown,  the  coverage  of  the  window  will  be  40  steps  in  the 
frequency  domain  and  64  steps  in  the  time  domain. 
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40xA/  -  2.5  64xAt  -  0.5  AREA  -   1.25 

where:         Af  =  0.0625  Hz      &      At  =  7.8125  msec 

The  AREA  =  1.25  is  a  unitless  quantity  that  defines  the  area  that  the  smoothing  window 
encompasses  while  smoothing  one  point  in  the  Distribution.  The  unitless  dimension  is  the 
result  of  the  time  being  in  seconds  and  the  frequency  in  hertz.  For  comparison  purposes,  the 
total  area  of  the  distribution  is  the  product  of  fmax  and  tmax  which  is  equal  to  256.0,  resulting 
in  the  window  smoothing  over  0.488%  of  the  total  Distribution,  while  smoothing  one  point. 
Now  consider  another  example  that  uses  1024  data  points.  As  shown  earlier,  the 
product  of  M  and  N  must  be  greater  than  or  equal  to  the  number  of  data  points  divided  by 
pi  in  order  to  meet  the  window  size  criterion  as  established  by  Reference  12.  With  respect  to 
the  earlier  example,  the  size  of  the  window  (ie.  N  times  M)  must  increase  two  times. 
Therefore,  a  selection  of  N  =  10  and  M  =  32  will  be  used  for  this  example.  These  values  are 
very  close  to  the  required  criterion  and  will  suffice  for  this  illustrative  example.  The  signal 
length  (tmax)  will  be  2.0  seconds  resulting  in  fmax  =  256.25  Hz,  At  =  1.953  msec,  and  Af  = 
0.125  Hz.   In  relation  to  the  area  of  coverage  of  the  window,  the  following  results  apply: 

40xAf-5.0  128xAt  -  0.249984  AREA  -  1.25 

where:         Af  =  0.125  Hz     &      At  =1.953  msec 

This  is  the  same  result  as  the  512  data  point  case,  with  the  following  exception.  The  total  area 
of  the  Distribution  (tmax  x  fmax  =  512.0)  has  doubled,  while  the  area  of  the  smoothing  window 
coverage  has  remained  the  same.  This  results  in  the  window  only  covering  0.244%  of  the 
entire  Distribution  when  smoothing  a  single  point  within  the  Distribution.  In  conclusion,  the 
greater  the  number  of  sampled  data  points  within  the  processed  raw  time  signal,  the  smaller 
the  area  the  smoothing  window  will  encompass  while  smoothing  the  Distribution.     An 
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important  trade-off  to  consider,  however,  is  the  increase  in  computer  time  required  to  run 
the  additional  number  of  data  points. 

D.       TIME/FREQUENCY  RESOLUTION 

Since  the  smoothing  window  used  in  this  work  gives  the  ability  to  vary  the  time  and 
frequency  size  of  the  window,  the  following  examples  show  the  results  that  can  be  obtained. 
Figure  5,  is  the  input  signal  that  is  sampled  using  512  data  points.  It  is  a  400  Hz  sine  wave 
that  has  been  computer  generated  with  an  amplitude  of  2.0.  In  order  to  be  consistent  with  the 
previous  work,  the  size  criterion  established  earlier  will  be  used  as  a  guide  for  the  selection 
of  M  and  N  in  the  smoothing  window.  Figure  6,  is  the  result  of  a  window  in  which  M  and 
N  are  both  equal  to  13.  Since  this  smoothing  window  does  not  require  M  and  N  to  be  equal, 
varying  these  parameters  yielded  some  intersting  results.  Figure  7,  shows  the  result  of  the 
input  signal  being  smoothed  with  a  35  by  5  window.  This  means  that  N  =  35  and  M  =  5,  or 
in  other  words,  the  window  is  longer  in  the  frequency  domain  than  it  is  in  the  time  domain. 
The  resulting  distribution  has  an  amplitude  that  is  much  more  constant  throughout  the  entire 
time  length,  but  the  amplitude  itself  has  been  reduced  considerably  as  compared  to  the  13  by 
13  example.  Figure  8  is  just  a  different  view  of  this  example,  showing  the  loss  of  frequency 
resolution  with  a  35  by  5  window.  Just  the  opposite  case  to  this  example  is  where  the  window 
used  for  smoothing  is  a  5  by  35  window.  In  this  case,  the  window  is  longer  in  the  time 
domain  than  the  frequency  domain  and  just  the  opposite  results  are  obtained.  Figure  9,  shows 
the  result  of  this  window,  while  Figure  10  gives  a  view  showing  the  degree  of  frequency 
resolution  that  can  be  obtained.  This  can  be  very  beneficial  in  the  monitoring  of  machinery, 
since  in  many  cases  the  ability  to  differentiate  between  signals  that  are  very  close  to  one 
another  is  paramount. 
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Figure  5.   400  Hz  Sine  Wave  Input 
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400  Hz  SINE  WAVE 

SMOOTHING  WINPOwTjgJCjia 

AMPLITUDE  =  2.0 


Figure  6.   Pseudo  Wigner-Ville  Distribution  (13  X  13  Window) 
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400  Hz  SINE  WAVE 

FREQUENCY  RESOLUTION  ;  35  X  5 

AMPLITUDE  =  2.0 
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Frequency  (Hz) 


1000.0 


Figure  8.   Frequency  Resolution  (35  X  5  Window) 
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400  Hz  SINE  WAVE 

SMOOTHING  WINDOW  ;  5  X  35 

AMPLITUDE  =  2.0 


Figure  9.   Pseudo  Wigner-Ville  Distribution  (5  X  35  Window) 
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Figure  10.   Frequency  Resolution  (5  X  35  Window) 
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E.        TRANSITIONAL  PEAK  PHENOMENON 

In  works  that  have  given  examples  of  the  Wigner-Ville  Distribution  for  transient 
signals,  the  point  at  which  the  signal's  frequency  is  changing  from  increasing  to  decreasing 
or  vice  versa  has  a  large  amplitude  increase  in  the  Wigner-Ville  Distribution.  For  example, 
Rossano  [Ref.  15]  shows  an  artificially  generated  linear  chirp  and  also  an  actual  pump 
signature  that  has  what  appears  to  be  a  large  concentration  of  energy  at  the  transition  point. 
In  order  to  use  the  Wigner-Ville  Distribution  for  transient  machinery  monitoring,  there  must 
be  an  explanation  or  way  to  quantify  this  "ghost  peak"  that  is  occuring  at  the  transition  point. 
To  date,  there  has  not  been  an  explanation  to  this  phenomenon,  however,  by  looking  at  the 
unsmoothed  Wigner-Ville  Distribution,  the  reason  for  this  is  easily  seen.  Figure  1 1  is  a  linear 
chirp  sine  wave  that  is  used  as  the  input  to  the  Wigner-Ville  Distribution.  Figure  12,  shows 
the  resulting  unsmoothed  Wigner-Ville  Distribution  of  this  linear  chirp  and  as  discussed  and 
shown  earlier,  the  cross-terms  can  be  seen  occuring  exactly  at  the  mid-point  of  the  conflicting 
frequencies  and  times.  These  cross-terms  build  until  they  reach  a  maximum  concentration 
at  the  transition  or  frequency  change  point.  Since  these  are  nothing  more  than  cross-terms, 
there  must  be  a  way  to  deemphasis  these  through  the  judicial  use  of  the  smoothing  window. 
Figures  13  through  16  show  the  result  of  varying  the  smoothing  window  from  a  point  of 
essentially  no  effect  to  entire  elimination  of  the  ghost  peak. 
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Figure  11.   Linear  Chirp  Sine  Wave  Input 


27 


LINEAR  CHIRP  SINE  WAVE 

UNSMOOTHED  DISTRIBUTION 

AMPLITUDE  =  1.0 


Tinr»e 


Figure  12.   Unsmoothed  Wigner-Ville  Distribution 
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LINEAR  CHIRP  SINE  WAVE 

SMOOTHED  :  5  X  35 

AMPLITUDE  =  1.0 


Figure  13.   Pseudo  Wigner-Ville  Distribution  (5  X  35  Window) 
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LINEAR  CHIRP  SINE  WAVE 

SMOOTHED  ;  13  X  13 

AMPLITUDE  ■  1.0 
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Figure  14.   Pseudo  Wigner-Ville  Distribution  (13  X  13  Window) 
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LINEAR  CHIRP  SINE  WAVE 

SMOOTHED  :  18  X  10 

AMPLITUDE  =  1.0 
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Figure  15.   Pseudo  Wigner-Ville  Distribution  (18  X  10  Window) 
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LINEAR  CHIRP  SINE  WAVE 

SMOOTHED  ;  35  X  5 

AMPLITUDE  =  1.0 


*) 


Figure  16.   Pseudo  Wigner-Ville  Distribution  (35  X  5  Window) 
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This  chapter  has  shown  the  benefits  that  may  be  gained  through  the  use  of  the 
smoothing  window.  The  frequency  resolution  and  amplitude  are  the  two  major  areas  affected 
by  the  window  changes.  The  last  area  presented,  showed  the  reason  for  the  large  amplitude 
increase  at  the  frequency  transition  point  and  how  these  "ghost  peaks"  may  be  smoothed  and 
to  what  extent  deemphasized  by  adjusting  the  smoothing  window  size.  Within  the  next 
chapter,  a  thorough  energy  analysis  of  the  discretized  Pseudo  Wigner-Ville  Distribution  will 
be  conducted. 
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IV.     WIGNER-VILLE  DISTRIBUTION  ENERGY  ANALYSIS 

A.        ENERGY  DEFINITION 

As  mentioned  earlier,  the  Wigner- Ville  Distribution  is  essentially  the  Fourier  transform 
of  the  auto-correlation  of  an  input  signal.  This  may  be  thought  of  as  the  spectral  density 
function  in  a  three  dimensional  sense.  The  total  energy  of  the  Wigner-Ville  Distribution  was 
listed  earlier  as: 


—  f  fwDF(u>,  t)  dodt  -  i|r 


where: 


1J/2  -  mean  square  value  of  the  signal 

There  are  a  few  points  that  must  be  made  concerning  this  equation.  The  first  is  that  the 
integration  is  from  minus  infinity  to  plus  infinty  and  another  is  that  it  is  the  integration  of 
the  resulting  Distribution  of  a  real  signal,  not  of  an  analytic  signal  as  used  in  most  discretized 
forms  of  the  Wigner-Ville  Distribution.  Lastly  is  that  this  equation  has  been  proven 
mathematically,  which  unfortunately  does  not  take  into  account  the  fact  that  there  are 
resulting  cross-terms  and  that  the  Distribution  is  subsequently  smoothed  because  of  this. 
Since  the  last  point  is  indicative  of  what  is  seen  in  actual  real  world  problems,  the  energy  must 
be  defined  in  a  discretized  sense  that  accounts  for  the  above  mentioned  points  and  also  these 
unavoidable  cross-terms. 
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B.        EFFECT  OF  WINDOW  SIZE  ON  ENERGY 

Since  smoothing  of  the  Distribution  is  a  necessary  requirement,  the  size  of  the  applied 
window  is  something  that  must  be  decided  upon.  Another  important  point  is  knowing  what 
effect  this  will  have  on  the  energy  of  the  Distribution.  In  the  previous  chapter,  the  size 
criterion  that  was  established  by  Cartwright  [Ref.  19]  and  later  used  by  Wahl  and  Bolton  in 
Reference  12  was  discussed.  Within  this  discussion,  the  visual  benefits  of  varying  the  window 
size  in  the  time  and  frequency  domain  were  presented.  In  order  to  define  the  energy,  the 
effect  on  the  energy  of  the  Pseudo  Wigner-Ville  Distribution  as  caused  by  varying  these 
windows  will  be  presented. 

Figure  17  is  a  400  Hz  sine  wave  of  length  0.256  seconds  and  amplitude  2.0  that  will  be 
sampled  and  have  the  Wigner-Ville  Distribution  calculated.  Figure  18  is  the  resulting 
smoothed  Wigner-Ville  Distribution  that  has  been  processed  using  a  13  by  13  sized  smoothing 
window.  Figure  18  also  shows  that  the  volume  is  equal  to  0.0709118  V2.  This  is 
representative  of  the  energy  contained  within  the  Wigner-Ville  Distribution  of  the  signal  in 
Figure  17  and  is  nothing  more  than  the  volume  under  the  resulting  Distribution.  By 
definition,  this  volume  should  be  equal  to  the  mean  square  value  of  the  400  Hz  sine  wave, 
however,  due  to  the  points  mentioned  above,  this  is  not  going  to  hold  true  and  will  be 
discussed  in  depth  later.  Figures  19  and  20  are  also  smoothed  Distributions  of  this  same  400 
Hz  signal  only  using  a  35  by  5  window  in  the  case  of  Figure  19,  and  a  5  by  35  window  in 
Figure  20.  Since  N  times  M  (the  size  criterion)  are  very  close  in  these  three  examples,  the 
resulting  energy  values  of  the  three  cases  are  essentially  the  same.  The  point  to  be  emphasized 
from  this  is  that  as  long  as  the  area  of  the  window  does  not  vary  significantly,  the  energy  will 
be  represented  in  a  constant  manner  and  independent  of  the  time  and  frequency  window  size 
selected.  The  last  example  that  will  be  presented  in  this  section  is  that  this  result  is 
independent  of  the  frequency  of  the  signal  as  shown  by  the  combined  100  and  400  Hz  sine 
wave  example  of  Figure  21,  in  which  the  energy  is  equal  to  0.1418314  V2.  This  is  essentially 
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twice  the  value  of  the  single  sine  waves  of  Figures  18  through  20.  Of  special  importance  for 
transient  machinery  applications  is  Figure  22,  which  is  a  transient  or  ramped  sine  wave  that 
shows  the  same  resulting  energy  as  the  previous  stationary  signals.  This  shows  that  the 
Wigner-Ville  Distribution  energy  is  independent  of  the  nature  of  the  signal  allowing  for  the 
possibility  of  monitoring  a  piece  of  machinery  during  both  it's  steady  state  and  transient 
phases  of  operation. 
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Figure  17.   400  Hz  Sine  Wave  Input 


37 


400  Hz  SINE  WAVE 

SMOOTHING  WINDOW  ;  13  X  13 

AMPLITUDE  =  2.0 


•  volume  ■  0.0709118  V~2 


Figure  18.   Pseudo  Wigner-Ville  Distribution  Energy  (13  X  13  Window) 
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400  Hz  SINE  WAVE 

SMOOTHING  WINDOW  J35  X  5 

AMPLITUDE  ■  2.0 


volume  -  0.0725799  \T2 


Figure  19.   Pseudo  Wigner-Ville  Distribution  Energy  (35  X  5  Window) 
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400  Hz  SINE  WAVE 

SMOOTHING  WINDOW  :  5  X  35 

AMPLITUDE  ■  2.0 


•  volume  -  0.0723785  V~2 


Figure  20.   Pseudo  Wigner-Ville  Distribution  Energy  (5  X  35  Window) 
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100  *  400  HZ  SINE  WAVE 

SMOOTHING  WINDOW  yJ3X  13 

AMPLITUDES  =  2.0 


•  volume  ■  0.1418314  \T2 


Figure  21.   Combined  Sine  Wave  Distribution  Energy  (13  X  13  Window) 
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RAMPED  SINE  WAVE 

SMOOTHING  WINDOW  ;  13  X  13 

AMPLITUDE  =  2.0 


•  volume  «  0.0709212  V~2 


Figure  22.    Ramped  Sine  Wave  Distribution  Energy  (13  X  13  Window) 
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The  previous  figures  show  that  if  the  Wigner-Ville  Distribution  is  going  to  be  used  to 
accurately  portray  or  maintain  signal  information,  then  the  window  size  must  at  least  remain 
constant  and  better  yet  to  keep  the  identical  window  size  for  each  measurement  of 
comparison.  This  requirement  to  maintain  a  constant  sized  smoothing  window  is  not  the  only 
constraint  that  must  be  meet  for  energy  analysis  purposes.  The  next  section  will  discuss 
another  of  these  elements. 

C.        TIME  DEPENDENCE  OF  WIGNER-VILLE  ENERGY 

Not  only  is  the  window  size  going  to  influence  the  energy  result,  but  the  time  length 
of  the  signal  will  also  play  a  large  part.  This  varies  considerably  from  the  conventional 
knowledge  of  the  mean  square  value  of  a  signal  being  independent  of  time.  The  standard 
mean  square  value  will  be  the  same  for  a  signal  of  time  length  2.0  seconds  as  it  will  for  a 
signal  of  4.0  seconds.   For  instance,  the  mean  square  value  of  a  sine  wave  x(t)  is: 


2 


t|r2  -   f(auto  spectrum)  df 
o 

where:         x(t)  ■  X  sin  (wt) 


To  compute  the  mean  square  value  of  this  sine  wave,  it  is  independent  of  time  and  will  equal 
0.5  for  X  =  1.0.  This  is  not  the  case  for  the  mean  square  value  or  total  energy  of  the  signal 
as  found  using  the  discretized  form  of  the  Wigner-Ville  Distribution.  The  following  figures 
will  show  this  second  important  point  that  must  be  realized  when  comparing  Wigner-Ville 
energies  of  various  signals.  Figure  23  is  the  400  Hz  input  signal  of  length  0.128  seconds.  This 
time  length  is  half  of  that  of  the  400  Hz  signal  used  in  the  previous  section.  The  resulting 
Wigner-Ville  Distribution  is  shown  in  Figure  24  with  an  energy  value  of  0.0354558  V2.  As 
can  easily  be  seen,  this  resulting  energy  value  is  not  the  same  and  in  fact,  the  energy  of  Figure 
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24,  is  one  half  that  of  the  earlier  400  Hz  signals  of  length  0.256  seconds.  Had  the  signal  been 
of  a  length  0.512  seconds  or  twice  as  long  as  the  earlier  examples,  the  resulting  Wigner-Ville 
energy  would  have  doubled.  The  window  size  used  for  this  Pseudo  Wigner-Ville  Distribution 
was  a  13  by  13  sized  window. 
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Figure  23.   400  Hz  Sine  Wave  (TL  =  0.128  sec) 
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400  HZ  SINE  WAVE 

SMOOTHING  WINDOW  ;  13  X  13 

AMPLITUDE  =  2.0 


•  volume  ■  0.0354558  \T2 


Figure  24.   Pseudo  Wigner-Ville  Distribution  Energy  (TL  =  0.128  sec) 
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This  has  yielded  a  second  concern  that  must  be  kept  in  mind  when  using  the  energy  of 
the  discretized  Pseudo  Wigner-Ville  Distribution.  Since  the  energy  is  proportional  to  the 
signal  length,  an  easy  conversion  can  be  made  by  using  the  lengths  of  the  individual  raw  input 
signals.  This  time  dependence  holds  true  for  not  only  steady  state  signals  as  shown,  but  also 
for  signals  of  a  transient  nature  such  as  previously  shown  with  Figure  22.  Throughout  these 
last  two  sections,  the  comparison  of  energy  values  of  signals  has  been  mentioned  a  great  deal. 
The  thought  process  behind  this  is  that  if  there  is  an  amplitude  change  in  the  input  signal,  this 
will  be  reflected  in  the  resulting  Pseudo  Wigner-Ville  Distribution.  By  ensuring  that  these 
previous  two  conditions  are  met,  an  investigation  into  the  resulting  energy  changes  from  an 
amplitude  increase  or  decrease  of  an  input  signal  can  be  made. 

D.        ENERGY  PROPORTIONALITY  OF  THE  DISTRIBUTION 

The  mean  square  value  of  a  simple  sine  wave  is  the  amplitude  squared  divided  by  two. 
This  was  shown  in  an  earlier  equation  and  now  will  be  used  to  show  the  proportionality  of  the 
Wigner-Ville  Distribution.  For  instance,  the  following  applies  for  signals  of  differing 
amplitudes  as  shown: 

amplitude  -  1.0  ijr2  -  0  . 5 

ampli  tude  -2.0  \|r2  -  2  .  0 

This  amplitude  increase  from  1.0  to  2.0  represents  an  increase  of  four  in  the  mean  square 
value  of  the  signal.  Figures  25  and  26  are  the  raw  input  sine  waves  of  amplitude  1.0  and  2.0 
respectively.  The  resulting  Pseudo  Wigner-Ville  Distributions  are  shown  in  Figures  27  and 
28,  which  acurately  portray  the  four-fold  energy  increase.  These  Distributions  have  been 
smoothed  using  a  window  of  13  by  13. 
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Figure  25.   400  Hz  Sine  Wave  (Amplitude  =  1.0) 
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Figure  26.   400  Hz  Sine  Wave  (Amplitude  =  2.0) 
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400  HZ  SINE  WAVE 

SMOOTHING  WINDOW  ;  13  X  13 

AMPLITUDE  =  1.0 


•  volume  -  0.0177279  \T2 


Figure  27.   Pseudo  Wigner-Ville  Distribution  (Amplitude  ■  1.0) 
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400  Hz  SINE  WAVE 

SMOOTHING  WINDOW  :  13  X  13 

AMPLITUDE  =  2.0 


•  volume  ■  0.0709118  \T2 


Figure  28.   Pseudo  Wigner-Ville  Distribution  (Amplitude  »  2.0) 
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These  previous  four  figures  show  the  proportionality  of  the  energy  within  the  Wigner- 
Ville  Distribution.  These  two  cases  are  just  an  example  to  show  the  proportionality  of  the 
energy.  Runs  were  made  for  many  more  cases  (ie.  X  =  4.0,  X  =  8.0,  etc)  and  the 
proportionality  relation  held  true  as  shown  with  the  previous  Figures.  This  was  another 
required  element  in  order  to  use  the  energy  as  the  basis  for  a  method  of  machinery 
monitoring.  Had  the  cross-terms  increased  in  such  a  fashion  as  to  not  allow  the  Distribution 
to  show  the  proportional  energy  change,  the  energy  content  of  the  Pseudo  Wigner-Ville 
Distribution  would  have  been  rendered  useless  for  the  purposes  of  comparison  and  evaluation. 
The  next  chapter  will  use  these  requirements  and  basics  that  have  been  developed  to  show  an 
actual  application  of  how  the  Wigner-Ville  energy  may  be  used  to  conduct  machinery 
monitoring. 
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V.     GEAR  MODEL  APPLICATION 

A.  DESCRIPTION  OF  GEAR  MODEL 

The  simple  complexity  gear  model  used  in  this  work  was  built  based  upon  a  model  used 
by  D.K.  Carlson  [Ref.  20]  for  Artificial  Neural  Network  applications  for  machinery 
monitoring.  It  may  be  of  a  simple  nature,  however,  it  provides  all  of  the  necessary 
components  that  may  be  found  in  more  complex  systems.  A  schematic  of  the  gear  model  can 
be  seen  in  Figure  29.  The  two  major  changes  of  this  model  over  the  prototype  used  by 
Carlson  is  that  the  components  are  much  larger  enabling  a  better  signal  to  noise  ratio  and 
secondly,  the  motor  used  to  drive  the  gear  model  is  mounted  on  a  separate  base  to  reduce 
motor  enduced  vibrations. 

The  model  consists  of  a  single  reduction  gear  train,  with  the  pinion  gear  being  a  90 
tooth  spur  gear  and  the  driven  gear  a  120  tooth  spur  gear.  The  gears  were  manufactured  by 
Martin  Corporation  with  a  1/2  inch  bore  and  a  14.5  degree  pressure  angle.  Aluminum  blocks 
housed  the  NICE  1/2  inch  bore  radial  ball  bearings  that  were  used  to  support  the  shaft  and 
gears.  These  aluminum  blocks  were  bolted  to  the  support  base  which  was  a  1.0  inch  thick 
plexiglass  slab.  Driving  the  single  shaft  is  a  General  Electric  l/6th  HP  motor,  that  is 
controlled  by  the  use  of  a  Bodine  Electric  Company  combination  rectifier  and  variable 
potentiometer  speed  controller.  The  shaft  speed  was  monitored  through  use  of  a  Power 
Instruments  Model  1720  RPM  Optical  Proximeter.  [Ref  20] 

B.  SIGNAL  SAMPLING  EQUIPMENT 

This  section  will  describe  the  components  shown  in  Figure  30,  that  were  used  to  sample 
the  vibrational  signal  of  the  gear  model.  The  signal  was  monitored  using  an  ENDEVCO 
model  303A03  accelerometer  that  was  amplified  and  powered  with  a  PCB  model  483B07 
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power  unit.  The  signal  was  then  processed  through  a  Krohn-Hite  model  3342  analog  filter 
configured  in  a  band  pass  mode,  prior  to  being  sent  to  the  HP  3565S  system  for  sampling  and 
storage.  Lastly,  the  signal  was  transfered  to  the  VAX  3520  workstation  for  the  calculation  and 
plotting  of  the  Pseudo  Wigner-Ville  Distribution. 

1.  ENDEVCO  Model  303 A03  Accelerometer 

The  ENDEVCO  Isotron  PE  accelerometer  is  a  minature  accelerometer  that  has  a 
medium  range  high  frequency  capability.  The  operation  of  the  accelerometer  is  based  upon 
a  piezoelectric  quartz  transducer  sensing  element.  The  sensitivity  of  the  accelerometer  is  10 
mV/g  with  a  resonant  frequency  of  70  kHz.  The  resolution  available  is  0.02  g,  with  a 
maximum  range  of  +  500  g. 

The  accelerometer  was  mounted  on  top  of  the  aluminum  bearing  housing  that 
supports  the  shaft  and  specifically  at  the  position  shown  in  Figure  30.  The  accelerometer  wass 
attached  by  means  of  a  mounting  wax  which  yields  a  maximum  operating  range  out  to  15k  - 
20k  Hz.   This  accelerometer  signal  was  then  sent  through  the  power  unit  described  below. 

2.  PCB  Model  483307  Power  Unit 

This  power  unit  is  used  to  power  the  low  impedance  quartz  transducers  and 
amplify  the  signal  if  desired.  The  power  unit  provides  an  adjustable  2  to  20  mA  current  for 
purposes  of  transducer  excitation.  The  gain  adjustment  is  available  from  0  to  100  and  is  set 
through  the  use  of  a  ten  turn  vernier  gain  pot  and  a  three  position  gain  multiplier  switch.  For 
the  purposes  of  these  tests,  the  gain  was  set  to  20  before  sending  the  signal  to  the  analog  filter. 

3.  Krohn-Hite  Model  3342  Anolog  Filter 

The  model  3342  variable  filter  is  a  digitally  tunned  filter  that  will  function  as  a 
High-Pass  or  Low-Pass  filter.  When  the  two  channels  of  the  filter  are  connected  together,  the 
filter  will  function  as  a  band  pass  filter,  which  is  the  configuration  for  the  gear  model  work. 
The  range  of  the  filter  is  from  0.001  Hz  to  99.9  kHz  as  adjusted  by  three  rotary  decade 


54 


switches  and  a  rotary  six  positiohn  multiplier  switch.    In  addition,  the  filter  unit  has  a  gain 
setting  of  unity  (0  dB)  and  10  (20  dB).    The  signal  was  not  further  amplified  using  this 
capability  of  the  filter  prior  to  sampkling  by  the  HP  3565S  system. 
4.       HP  3565S  Measurement  Hardware  System 

The  HP  3565S  measurement  hardware  with  the  HP- Vista  software  uses  a  HP-9000 
series  computer  system  for  controlling  purposes.  The  HP  3565S  system  is  a  modular 
multichannel  system  that  can  process  data  in  both  the  time  and  frequency  domain.  The 
system  is  capable  of  handling  64  source  and  input  modules,  but  is  presently  configured  for 
eight.  The  HP- Vista  software  allowed  for  the  sampling  and  storage  of  the  preprocessed 
accelerometer  signal  that  was  later  processed  through  the  Pseudo  Wigner-Ville 
distribution. 
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Figure  29.  Gear  Model  Schematic 
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Figure  30.  Signal  Sampling  Equipment  Schematic 
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C.        ESTABLISHING  MONITORING  FREQUENCIES 

This  section  will  describe  how  the  monitoring  frequencies  have  been  determined  and 
how  they  will  be  used  to  monitor  the  gear  model  and  in  future  applications  a  piece  of 
transient  machinery.  As  shown  earlier  in  Chapter  III,  the  energy  contained  in  a  signal  is 
accurately  portrayed  within  the  resulting  time-frequency  domain  of  the  Wigner-Ville 
Distribution.  Through  the  Wigner-Ville  Distribution  energy  that  is  contained  in  the 
established  monitoring  frequencies,  faults  may  be  detected.  The  shaft  frequency  is  just  the 
speed  of  the  shaft,  while  the  gear  mesh  frequency  is  the  number  of  teeth  on  the  rotating  gear 
multiplied  by  the  shaft  speed  or  frequency.  Both  of  these  frequencies  have  associated 
harmonics  that  can  be  of  equal  importance  as  the  frequencies  themselves  in  the  monitoring 
of  machinery,  It  is  for  these  reasons  that  the  frequencies  listed  below  have  been  chosen  based 
upon  a  shaft  frequency  of  15  Hz. 

1.       Monitoring  Frequencies 

For  the  purposes  of  this  work,  just  the  shaft  and  gear  mesh  frequencies  of  the 
model  will  be  included.  Just  as  in  standard  established  machinery  monitoring  methods,  the 
shaft  frequency  and  associated  harmonics  can  provide  a  great  deal  of  information  for  fault 
detection  and  diagnosis.  Along  with  the  shaft  frequencies,  the  gear  meshing  frequencies 
provide  valuable  information  of  faults  and  will  also  be  monitored.  The  purpose  of  the  two 
examples  will  be  to  show  that  the  Wigner-Ville  Distribution  energy  (ie.  volume)  changes  as 
a  result  of  damage  and  will  prove  to  be  a  very  valuable  tool  in  future  machinery  monitoring 
methods.  The  established  monitoring  frequencies  have  been  filtered  as  follows: 

•  Shaft  Input  1:  (5-100  Hz)  This  is  the  spectrum  of  signals  in  the  lower  frequency  range 
selected  to  determine  the  total  energy  contained  within  the  shaft  frequency  and  it's 
harmonics. 

•  Shaft  Input  2:  (10-20  Hz)  This  frequency  range  captures  just  the  shaft  frequency. 

•  Shaft  Input  3:  (25  -  35  Hz)  This  frequency  range  captures  the  l8t  shaft  harmonic. 
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•  Shaft  Input  4:   (40  -  50  Hz)   This  frequency  range  captures  the  2nd  shaft  harmonic. 

•  Shaft  Input  5:  (55  -  100  Hz)  This  frequency  range  captures  the  upper  harmonics  of  the 
shaft  frequency. 

•  Gear  Mesh  Input  1:  (1250  -  1450  Hz)  This  frequency  range  captures  the  gear  meshing 
frequency. 

•  Gear  Mesh  Input  2:  (2600  -  2800  Hz)  This  frequency  range  captures  the  1st  gear  mesh 
harmonic. 

•  Gear  Mesh  Input  3:  (3950  -  4150  Hz)  This  frequency  range  captures  the  2nd  gear  mesh 
harmonic. 


2.       Combined  Monitoring  Technique 

Using  the  Wigner-Ville  energy  as  an  input  to  an  Artificial  Neural  Network  is  a 
follow-on  goal  of  this  research.  It  is  for  this  reason,  that  the  above  established  frequencies 
have  been  identified  as  inputs.  Within  this  scheme,  the  Wigner-Ville  is  used  to  characterize 
the  signal  and  give  a  single  energy  value  at  the  particular  frequency  of  interest  as  an  input  to 
an  Artificial  Neural  Network.  By  filtering  these  signals  prior  to  processing,  just  the 
frequency  of  interest  will  be  captured  along  with  it's  respective  Wigner-Ville  energy  content. 
The  Wigner-Ville  energy  data  can  now  be  used  to  produce  a  comprehensive  training  set  for 
the  Artificial  Neural  Network.  By  pre-processing  the  inputs  to  a  Neural  Network  with  the 
filtering  and  the  Wigner-Ville  Distribution,  the  number  of  inputs  is  greatly  reduced.  Since 
these  associated  energy  levels  are  fairly  stable  for  a  given  machinery  condition,  the  Neural 
Network  will  be  capable  of  detecting  patterns  in  energy  level  shifts  that  are  associated  with 
particular  machinery  faults.  To  support  this  last  statement,  the  following  section  will  show 
two  examples  of  machinery  faults  imposed  on  the  gear  model  and  the  response  of  the  Wigner- 
Ville  energy  to  these  faults. 
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D.        MACHINERY  FAULT  EXAMPLES 
1.       Mechanical  Looseness 

This  fault  was  imposed  on  the  gear  model  by  loosening  the  aluminum  support 
blocks  that  house  the  bearings  and  support  the  shaft.  Specifically,  the  blocks  that  support  the 
driver  shaft  were  loosened.  The  degree  of  looseness  imposed  was  difficult  to  judge  and  made 
this  particular  fault  a  difficult  one  to  diagnose.  The  blocks  were  not  completely  undone,  but 
loosened  just  a  small  amount  from  their  original  firmly  bolted  posture.  Shown  in  the  table 
below  are  the  results  of  averaging  two  good  runs  and  two  post-damage  runs  and  comparing 
the  respective  energy  values  at  the  different  frequencies.  These  energy  values  are  nothing 
more  than  the  volume  under  the  Pseudo  Wigner-Ville  Distribution  as  discussed  earlier.  Only 
the  shaft  frequencies  are  shown  for  purposes  of  clarity  and  since  it  is  in  these  monitored 
frequencies  that  a  fault  of  mechanical  looseness  should  present  itself. 
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Table  I.   MECHANICAL  LOOSENESS 


MECHANICAL  LOOSENESS 

ENERGY  VALUES  (V~2) 


Freauency 

Good 

Damaaed 

Shaft  Spectrum 

41.0 

28.0 

Shaft  Frequency 

1.07 

0.83 

1st  Harmonic 

2.31 

1.22 

2nd  Harmonic 

1.23 

2.40 

Upper  Spectrum 

31.4 

26.31 

As  can  be  seen  in  the  table  above,  the  energy  content  approximately  doubled  at  the  second 
shaft  harmonic.    This  is  a  very  encouraging  result  for  the  purposes  of  using  the  Pseudo 
Wigner-Ville  Distribution  energy  as  a  fault  detection  tool. 
2.      Gear  Tooth  Damage 

This  particular  fault  was  imposed  on  the  gear  model  by  shaving  a  gear  tooth  on 
the  driver  gear.  The  degree  of  damage  was  much  easier  to  predict  in  this  case  and  made  this 
particular  fault  an  easier  one  to  control.  The  gear  tooth  that  was  damaged  had  approximately, 
l/8th  of  the  leading  edge  removed.  Shown  in  the  table  below  are  again  the  results  of 
averaging  two  good  runs  and  two  post-damage  runs  and  comparing  the  respective  energy 
values  at  the  different  frequencies.  For  this  case,  just  the  frequencies  associated  with  the  gear 
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mesh  are  presented,  since  it  is  in  these  monitored  frequencies  that  a  fault  involving  the 
damage  of  a  gear  train  would  present  itself. 
Table  II.   GEAR  TOOTH  DAMAGE 

GEAR  TOOTH  DAMAGE 

ENERGY  VALUES  (V~2) 

MESH  FREQUENCY      GOOD     DAMAGED 


Fundamental 

17.14 

154.4 

1st  Harmonic 

0.53 

3.44 

2nd  Harmonic 

0.16 

0.48 

As  shown  in  the  above  table,  the  energy  content  greatly  increased  at  the  fundamental  gear 
mesh  frequency.  Additionaly,  there  was  also  increases  at  the  harmonics  and  shown  in  Figures 
31  and  32  is  a  graphical  presentation  of  the  energy  increase  that  occured  due  to  the  damaged 
tooth. 
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GEAR  TOOTH  DAMAGE 
GEAR  MESH  INPUT  ONE 
FILTERED  1250  -  1450  HZ 


** 


Figure  31.  Gear  Mesh  Frequency  (Pre-Damage) 
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GEAR  TOOTH  DAMAGE 
GEAR  MESH  INPUT  ONE 
FILTERED  1250  -  1450  HZ 


Figure  32.   Gear  Mesh  Frequency  (Post-Damage) 
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This  chapter  has  shown  two  examples  of  machinery  damage  and  the  resulting  changes 
of  the  Pseudo  Wigner-Ville  Distribution  energy.  By  applying  the  window  and  energy 
requirements  discussed  in  Chapters  III  and  IV,  the  Wigner-Ville  energy  has  yielded  a  definite 
change  to  allow  a  means  for  machinery  monitoring.  The  Wigner-Ville  energy  will  make  a 
stable  input  to  a  Neural  Network  and  provide  a  quick  efficient  method  by  pre-processing  the 
machinery  signal  with  the  Wigner-Ville  Distribution  as  shown.  This  pre-processing  will  allow 
the  number  of  inputs  to  the  Neural  Network  to  remain  at  a  reasonable  level  and  provide  for 
a  manageable,  quicker  hardware  system. 
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VI.     CONCLUSIONS 

The  Pseudo  Wigner-Ville  Distribution  has  been  shown  to  provide  an  excellent 
characterization  of  an  input  signal.  Based  upon  the  research  conducted,  it  will  also  provide 
a  means  for  machinery  condition  monitoring  and  diagnostics.  The  research  has  also  shown 
that  in  order  to  use  the  Wigner-Ville  Distribution  in  a  machinery  monitoring  scheme,  the 
following  conditions  or  conclusions  must  be  met  and  understood. 

•  The  smoothing  window  size  applied  must  remain  constant  for  purposes  of  comparing 
various  signals. 

•  The  individual  time  and  frequency  lengths  of  the  smoothing  window  may  be  adjusted 
to  provide  a  better  frequency  resolution  or  constant  amplitude. 

•  The  energy  contained  within  the  Pseudo  Wigner-Ville  Distribution  of  a  signal  is 
proportionaly  time  dependent. 

•  The  energy  of  the  Distribution  is  proportional  to  the  mean  square  value  of  the  input 
signal  and  is  not  adversely  affected  by  the  cross-terms. 

•  Lastly,  the  damaged  gear  model  conditions  showed  a  definite  change  in  the  energy  of 
the  Pseudo  Wigner-Ville  Distribution  at  the  respective  monitoring  frequencies. 
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VII.     RECOMMENDATIONS  FOR  FUTURE  RESEARCH 

The  Pseudo  Wigner-Ville  Distribution  computer  program  that  was  established  here  at 
the  Naval  Postgraduate  School  is  very  stable.  The  pronounced  peaks  at  the  transition  points 
are  now  definable  and  additionally,  with  the  energy  of  the  Pseudo  Wigner-Ville  Distribution 
now  being  a  definable  and  understood  quantity,  the  next  additional  research  in  this  area  may 
include: 

•  Establish  a  baseline  or  data  bank  of  "good"  Wigner-Ville  Distribution  data  for  the 
present  gear  model. 

•  Conduct  additional  runs  of  various  levels  of  damage  of  gear  model  components  for 
comparison  with  good  data. 

•  Establish  a  Torpedo  Ejection  Pump  baseline  of  good  data  as  with  the  gear  model. 

•  Continue  with  the  research  of  linking  the  Pseudo  Wigner-Ville  Distribution  as  an  input 
to  a  Neural  Network  Application. 


66 


APPENDIX      COMPUTER  PROGRAMS 

The  following  programs  have  all  been  developed  here  at  the  Naval  Postgraduate  School 
and  used  in  the  conduct  of  the  research  presented.  The  language  used  is  FORTRAN  77  with 
the  plotting  programs  using  calls  to  CA-DISSPLA  software. 

A.   WVD3  COMPUTER  PROGRAM 

PROGRAM  WVD3 

* 
************************************************************ 

*  * 

*  THE  FOLLOWING  INDIVIDUALS  HAVE  BEEN  INSTRUMENTAL  IN  THE  * 

*  DEVELOPMENT  OF  THIS  PROGRAM  AND  PREVIOUS  VERSIONS:  * 

*  * 

*  * 

*  Graham  Rossano  * 

*  * 

*  Jose  G.  Calahorrano  * 

*  * 

*  Scott  G.  Spooner  * 

*  * 

*  Prof.  Y.S.  Shin  * 

*  * 

*  * 
************************************************************ 
* 

*  This  program  calculates  the  Wigner  Distribution 

*  Function  of  a  signal.   It  includes  the  option  of  applying  a 

*  smoothing  function  and  reducing  the  distribution  to  different 

*  sizes  while  writing  the  results  to  files  for  subsequent 

*  plotting. 

* 
************************************************************ 

*  * 

*  VERSION  3  * 

*  * 
************************************************************ 

*  • 

*  VARIABLE  DECLARATION  * 

*  * 
************************************************************ 
* 

**********...  setting  number  of  data  points  ---************* 
* 

* 

integer  dp,dp2,mvopt,redopt,dimt,dimf ,mm,nn,rf ,rdp,rdp2, 
:       nf ,mt,nf2,mt2 

PARAMETER   (DP  =  512) 

real  pi ,tin(dp),ain(dp),rawa(dp),rwdf (dp*2,dp),rswdf (dp*2,dp), 
:     rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del1 ,del2, 
:     const, t, sum, sumb,val,sval,wdf (dp*2,dp),coef ,df .AREA 

complex  s(dp*2),dum,c(dp*5),dum3,dum2 

character*25     inname 
* 

common/ in/  mvopt,redopt,dimt,dimf ,mm,nn,rf ,rdp,rdp2, 
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:       nf ,mt,nf2,mt2 
common/rl/  pi,tin,ain,rawa,rwdf , rswdf , 
:     raut ,dtsum,del t,dt ,asum,meanv,mt ime,del1 ,del2, 
:     const , t , sum, sumb, va I , sva I , wdf , coef , df , area 
common/ comp/  s,dum,c,dum3,dum2 
dp2  =  dp*2 
* 
************************************************************ 

*  INITIALIZING  VARIABLES  * 

************************************************************ 

*  —  Print  description  of  program  — 
print* 

print*,'  Program  to  calculate  the' 
print*,'  Pseudo  Wigner-Ville  Distribution' 
print* 

*  ---  Set  initial  values  --- 

print*,'  Enter  name  of  signal  input  file' 

read(5,904)  inname 

print* 

* 

print*,'  Do  you  wish  to  remove  the  mean  value?' 
print*,'  Enter  1  for  yes  -  0  for  no' 
read(5,910)  mvopt 
print* 


print*, ' 

Input  the  reduction  size  desired 

print*, ' 

input  1  for 

64  by  32       ' 

print*, ' 

input  2  for 

128  by  64      • 

print*, ' 

input  3  for 

128  by  128     ' 

print*, ' 

input  4  for 

256  by  128     ' 

read(5,910)  redopt 
print* 

* 

print*,'  Input  the  smoothing  window  size  desired' 

print*,'  input  the  size  for  the  frequency  axis' 

read(5,910)  nf 

print*,'  input  the  size  for  the  time  axis' 

read(5,910)  mt 
* 
* 
* 
a*********************************************************** 

*  * 

*  calculation  of  some  common  used  constants         * 

*  * 
************************************************************ 

pi  =  4.0  *  atan(LO) 


**************************************************************** 

*      THIS  IS  THE  CALCULATION  PART  OF  THE  PROGRAM  * 

**************************************************************** 


open(4,f i le=inname,status='old' ) 
rewind  4 


*  PRELIMINARY  CALCULATIONS: 
CALL  DATAIN3 
write(6,*)  ,  '  done  with  datain3' 
CALL  DTCALC3 
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write(6,*)  ,  '  done  with  dtcalc3' 

CALL  ME AN 3 
write(6,*)  ,  '  done  with  mean3' 
write(6,*)  ,  '  mean  value  removed  =',  meanv 

*  SIGNAL  MODIFICATIONS: 

*  Window  application  (modified  hamming  window) 

CALL  HAMMG3 
write(6,*)  ,  '  done  with  hammg3' 

*  Conversion  of  real  signal  to  an  analytical  one 

CALL  ANAL3 
write(6,*>  ,  '  done  with  anal3' 


*  CALCULATION  OF  THE  UIGNER  VILLE  DISTRIBUTION: 

CALL  UIGH3 
write(6,*)  ,  '  done  with  wigh3' 

*  REDUCTION  &  SMOOTHING  OF  THE  UIGNER  VILLE  DISTRIBUTION: 

CALL  REDUCE3 

write(6,*)  ,  ■  done  with  reduce3' 
CALL  SMOOTH3 

write(6,*)  ,  '  done  with  smooth3' 

*********************************************************** 

*  * 

*  WRITING  OF  THE  DIFFERENT  ARRAYS  GENERATED        * 

*  TO  FILES  FOR  PLOTTING  PURPOSES.  * 

*  * 
*********************************************************** 
* 

*  OPENING  OF  THE  OUTPUT  FILES 

* 

OPEN(UNIT=7,FILE=,RAWTIME.OUT'fSTATUS=,NEW') 

OPEN(UNIT=8,FILE=,WNDWTIM.OUT,fSTATUS=,NEW') 

OPEN(UNIT=9,FILE='RWDF.OUT,,STATUS='NEWl) 

OPEN(UNIT=10JFILE=,RSWDF.OUT'fSTATUS='NEW) 

OPEN(UNIT=11,FILE=,WDF.OUT',STATUS='NEWl) 


*  WRITING  OF  RAW  AND  WINDOWED  TIME  SIGNAL  TO  FILES 

* 

DO  300  I  =  1,DP 

WRITE(7,930)  TIN(I),RAWA(I) 
WRITE(8,930)   T!N( I  ),AIN( I  ) 

300     CONTINUE 

* 

*  WRITING  OF  REDUCED/UNSMOOTHED  WVD  TO  FILE 

* 

DO  400  I  =  1,DP/mm 
DO  400  J  =  1,DP2/nn 

WR!TE(9,906)   RWDF(J,I) 

400       CONTINUE 

* 

*  WRITING  OF  REDUCED  &  SMOOTHED  WVD  TO  FILE 
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DO  500  I  =  1,DP/mm 
DO  500  J  =  1,DP2/nn 

URITE<10,906)  RSWDF(J,I> 
500       CONTINUE 

*  WRITING  OF  ENTIRE  UDF  TO  A  FILE 

* 

DO  600  I  =  1.DP 
DO  600  J  =  1.DP2 

WRITE(11,906)   WDF(J.I) 
600       CONTINUE 

*  Format  statements 

904  format (a25) 

906  format(2X,g16.8) 

910  format<i"6) 

930  format(2x,g16.8,5x,g16.8) 

END 

*  SUBROUTINES 


include 
include 
include 
include 
i  nc  I  ude 
include 
include 
include 
include 
include 


•ANAL3.INC1 

•C0RR3.INC 

•DATAIN3.INC 

'DTCALC3.INC 

'FFT3.INC 

'HAMMG3.INC' 

'HEAN3.INC 

'SM00TH3.INC' 

'REDUCE3.INC 

■U1GH3.INC 
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B.        SUBROUTINES  OF  WVD3 


Listed  below  in  alphabetic  order  are  the  subroutines  used  in  the  main  program  WVD3. 


SUBROUTINE  ANAL3 


*  This  subroutine  converts  a  real  signal  to  an         * 

*  analytic  one  * 

*  * 
********************************************************** 


integer  dp,dp2,mvoptfredoptfdimt,dimf ,mm,nn,rf ,rdp,rdp2, 
:       nf ,mt,nf2,mt2 

PARAMETER   (DP  =  512) 

real  pi,tin(dp),ain(dp),rawa(dp),rudf (dp*2,dp),rswdf (dp*2,dp), 
:     rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del1 ,del2, 
:     const , t , sum, sumb, va I , sva I , wdf (dp*2 , dp) , coef , df .area 

complex  s(dp*2),dum,c(dp*5),dum3,dum2 

character*25  inname 

common/ in/  mvopt,redopt,dimt,dimf ,mm,nn,rf , rdp,rdp2, 
:       nf ,mt,nf2,mt2 

common/rl/  pi , tin,ain,rawa,rwdf ,rswdf , 
:     rawt, dt sum, del t,dt,asum,meanv, mtime.de 11 ,del2, 
:     const , t , sum, sumb, va I , sva I , wdf , coef , df , area 

common/comp/  s,dum,c,dum3,dum2 

dp2  =  dp*2 


do  100  i=1,dp 
sum=0.0 

do  200  j  =  1 , dp 
sumb=0.0 

if(i-j.eq.O)  go  to  200 

n=i-j 

val=pi*n/2. 

sval=sin(val) 

sumb=ain( j )*sval*sval/val 

200      sum=sum+sumb 

s(i)=cmplx(ain(i),sum) 

100      continue 

END 
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SUBROUTINE  DATAIN3 

************************************************************** 

*  * 

*  This  subroutine  assumes  the  input  file  is  in  the  following  * 

*  format:  time  amplitude  (2x,e16.8,5x,e16.8)  with  an  end  of  * 

*  file  indicator  of  a  last  entry  9999., 9999.  * 

*  * 
************************************************************** 

* 
* 
* 

integer  dp,dp2,mvopt,redopt,dimt,dimf ,mm,nn,rf ,rdp,rdp2, 
:       nf ,mt,nf2,mt2 

PARAMETER  (DP  =  512) 

real  pi , tin(dp),ain(dp), rawa(dp), rwdf ( dp*2, dp) , rswdf (dp*2,dp), 
:     rawt (dp), dtsum, del tfdt,asum,meanv,mtime, deM.de 12, 
:     const , t , sum, sumb, va I , sva I , wdf (dp*2 , dp) , coef , df , ARE A 

complex  s(dp*2),dum,c(dp*5),dum3,dum2 

character*25  inname 
* 

common/ in/  mvopt,redopt,dimt,dimf ,mm,nn,rf ,rdp,rdp2, 
:       nf ,mt,nf2,mt2 

common/r I /  pi,tin,ain,rawa,rwdf , rswdf , 
:     rawt, dt sum, del t,dt,asum,meanv,mt ime.de 11, del  2, 
:     const , t , sum, sumb, va I , sva I , wdf , coef , df , area 

common/comp/  s,dum,c,dum3,dum2 

dp2  =  dp*2 


********simple  loop  to  read  in  time  &  amplitude************* 

do  100  j  =  1 , dp 

read(4,902)  rawt(j)  ,  rawa(j) 
tin(j)  =  rawt(j) 
ain(j)  =  rawa(j) 

100     continue 

*    FORMAT  STATEMENT 

902     FORMAT(2X,E16.8,5X,E16.8) 

200      END 
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SUBROUTINE  DTCALC3 


*  This  subroutine  calculates  the  delta  t  of  the  signal 


integer  dp,dp2,mvopt,redopt,dimt,dimf ,nm,nn,rf ,rdp, rdp2, 
:       nf ,mt,nf2,mt2 
PARAMETER   (DP  =  512) 
real  pi , tin(dp),ain(dp),rawa(dp),rwdf (dp*2,dp),rswdf(dp*2,dp), 

:     rawt(dp),dtSLfn,del  t  ,dt  ,asLfti,meanv,mtime,del  1  ,del2, 
:     const , t , sum, sumb, va I , sva I , wdf (dp*2 , dp ) , coef , df , AREA 

complex  s(dp*2),dum,c(dp*5),dum3,dum2 

character*25  inname 

common/ i n/  mvopt , redopt , dimt ,di  mf , mm, nn, r f , rdp, rdp2, 
:       nf ,mt,nf2,mt2 

common/rl/  pi ,tin,ain,rawa,rwdf ,rswdf , 
:     rawt,dtsum,delt,dt,asum,meanv,mtime,del1 ,del2, 
:     const , t , sum, sumb, va I , sva I , wdf , coef , df , area 

common/comp/  s,dum,c,dum3,dum2 

dp2  =  dp*2 


* 
************************************************************ 

dtsum  =  0.0 

do  100  i  =  1  ,  dp-1 

delt  =  tin(i+1)  -  tin(i) 
dtsum  =  dtsum  +  delt 

1 00    cont  i  nue 

dt  =  dtsum  /  (dp-1.) 


END 
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SUBROUTINE    FFT3 

* 
********************************************************** 

*  * 

*  This  subroutine  calculates  the  Fast  Fourier  Transform  * 

*  * 
********************************************************** 


integer  dp,dp2,mvopt,redopt,dimt,dimf ,mm,m,rf , rdp,rdp2, 
:       nf ,mt,nf2,mt2 

PARAMETER   (DP  =  512) 

real   pi,tin(dp),ain(dp), rawa(dp), rwdf (dp*2,dp),rswdf (dp*2,dp), 
:  rawt(dp),dtsun,del t,dt , asum,meanv,mt ime.det 1 , del2, 

:  const , t , sum, suinb, va I , sva I , wdf (dp*2, dp) , coef , df 

complex  s(dp*2),dum,c(dp*5),dum3,dum2 

character*25  irmame 

common/ in/  mvopt,redopt,dimt,dimf ,mm,nn,rf , rdp,rdp2, 
:       nf ,mt,nf2,mt2 

common/rl/  pi,tin,ain,rawa,rwdf ,rswdf , 
:     rawt , dtsum.de I t,dt,  asum.meanv.mtime.dell.de^, 
:     const , t , sum, sumb, va I , sva I , wdf , coef , df , area 

common/ comp/  s , dum ,  c , duri5 , dun2 

dp2  =  dp*2.0 


const=dp2 

val=alog(const)/alog(2.)+.1 

j=1 

do  40  i=1,dp2-1 

if  (i.ge.j)  go  to  10 

dum3=c( j) 

c(j)=c(i) 

c(i )=dum3 
10     k=dp 
20     if  (k.ge.j)  go  to  30 

j=j-k 

k=k/2 

go  to  20 
30     j=j+k 
40   cont  i  nue 

do  70  n=1,val 

coef=2**n 

coef=coef/2 

dum2=cmplx(1.>0.) 

dum=cmpl x( cos(pi /coef ),- s i n(pi /coef ) ) 

do  60  j=1,coef 

do  50  i=j,dp2,2*coef 
ii=i+coef 
dum3=c(i  i )*dum2 
c(ii)=c(i)-dum3 
c(i)=c(i)+dum3 
50      continue 

dum2=dum2*dum 
60    cont  i  nue 
70   cont  i  nue 


END 
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SUBROUTINE  HAHHG3 
* 
******************************************************** 

*  * 

*  This  subroutine  applies  a  modified  hamming  window   * 

*  to  the  signal  ain(t)  * 

*  * 
******************************************************** 
* 

* 

i nt eger  dp, dp2 , mvopt , redopt , di mt , d i mf , mm, nn, r f , rdp, rdp2 , 
:       nf ,mt,nf2,mt2 

PARAMETER   (DP  =  512) 

real  pi , tin(dp),ain(dp),rawa(dp),rwdf (dp*2,dp),rswdf (dp*2,dp), 
:     rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
:     const , t , sum, sumb, va I , sva I , wdf (dp*2, dp) , coef , df , area 

complex  s(dp*2),dum,c(dp*5),dum3,dum2 

character*25  irmame 
* 

common/ i  n/  mvopt , redopt , d  i  mt , di  mf , mm, nn, rf , rdp, rdp2 , 
:       nf ,mt,nf2,mt2 

common/rl/  pi,tin,ain,rawa,rwdf ,rswdf , 
:     rawt , dt sum, de I t , dt , asum, meanv, mt  i  me , de 1 1 , de 1 2 , 
:     const , t , sum, sumb, va I , sva I , wdf , coef , df , area 

common/comp/  s,dum,c,dum3,dum2 

dp2  =  dp*2 


mtime=dp*dt 
del1=0.1*mtime 
del2=0.9*mtime 
const=pi/del1 

do  100  j  =  1  ,  dp 

t  =  (j-1)  *  dt 

if  (t.le.deM)  then 

ain(j)  =  ain(j)  *  ( .5*0 .-cos(const*t))) 

elseif  ((t.ge.del2).and.(t.le.mtime))  then 

ain(j)=ain( j)*(.5*(1.-cos(const*(mtime-t)))) 

else 

end  if 

100     continue 

END 
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SUBROUTINE   ME AN 3 
* 
********************************************************** 

*  * 

*  This  subroutine  calculates  and  removes  the  mean  value  * 

*  of  the  signal.  * 

*  * 
********************************************************** 


integer  dp,dp2,mvopt,redopt,dimt,dimf ,inm,nn,rf , rdp,rdp2, 
:       nf ,mt,nf2,mt2 

PARAMETER  (DP  =  512) 

real  pi\tin(clp),ain(dp),rawa(dp),rwdf(dp*2,dp),rswdf(dp*2,dp), 
:     rawt(dp),dtsun,del t ,dt, asum, meanv, mt ime.de  I 1 .del  2, 
:     const , t , sum, sumb, va I , sva I , wdf (dp*2 , dp) , coef , df , AREA 

complex  s(dp*2),dum,c(dp*5),dum3,dum2 

character*25  imame 

common/ in/  mvopt,redopt,dimt,dimf ,mm,nn,rf ,rdp,rdp2, 
:       nf ,mt,nf2,mt2 

common/rl/  pi ,tin,ain,rawa,rwdf .rswdf , 
:     rawt , dt  sum, de I t , dt , asum, meanv , mt  i  me , de 1 1 , de 1 2 , 
:     const , t , sum, sumb, va I , sva I , wdf , coef , df , area 

common/comp/  s,dum,c,dum3,dum2 

dp2  =  dp*2 


asum  =0.0 

do  100  i  =  1  ,  dp 

asum  =  asum  +  ain(i) 

100     continue 

meanv  =  asum  /  dp 

if  (mvopt.eq.1)  then 

do  200  i  =  1  ,  dp 
ain(i)  =  ain(i)  -  meanv 
200        continue 

endif 

END 
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SUBROUTINE  REDUCE3 


* 

*  This  subroutine  reduces  the  wigner  ville 

*  to  a  workable  size. 

* 


i  nteger  dp, dp2 , mvopt , redopt , di mt , d i mf , mm, nn, rf , rdp, rdp2 , 
:       nf ,mt,nf2,mt2 

PARAMETER   (DP  =  512) 

real  pi ,tin(dp),ain(dp),rawa(dp),rwdf (dp*2,dp),rswdf (dp*2,dp), 
:     rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
:     const , t , sun, sumb, va I , sva I , wdf (dp*2, dp) , coef ,df , area 

complex  s(dp*2),dum,c(dp*5),dum3,dum2 

character*25  irmame 

common/ in/  mvopt, redopt, dimt.dimf , mm, rm,rf , rdp, rdp2, 
:       nf ,mt,nf2,mt2 

common/ r I /  pi , t  i  n, a i  n , rawa , r wdf , rswdf , 
:     rawt , dt  sum, de 1 1 , dt , asum, meanv, mt  i  me , de 1 1 , de 1 2 , 
:     const , t , sum, sumb, va I , sva I , wdf , coef , df , area 

common/comp/  s,dum,c,dum3,dum2 

dp2  =  dp*2 


IF  (REDOPT. EQ.1)  THEN 
RF  =  64 
RT  =  32 
ELSEIF  (REDOPT. EQ. 2)  THEN 
RF  =  128 
RT  =  64 
ELSEIF  (REDOPT. EQ. 3)  THEN 
RF  =  128 
RT  =  128 
ELSE 

RF  =  256 
RT  =  128 
END  IF 

nn  =  dp2  /  RF 
mm  =  dp  /  RT 

I  =  0 
do  100  j  =  1  ,  dp2  ,  nn 
1  =  1  +  1 
k  =  0 
do  100  i  =  1  ,  dp  ,  mm 
k  =  k  +  1 

rwdf(l.k)  =  wdf(j.i) 
100    continue 

END 
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SUBROUTINE  SMOOTH3 


*  This  subroutine  reduces  and  smooths  the  UDF  along  both    * 

*  the  time  and  frequency  axes  for  purposes  of  plotting     * 

*  clarity.  * 

*  * 

************************************************************* 


* 


* 


* 


real  hg(-100:100, -100:100) 


integer  dp.dpZ.mvopt.redopt.dimt.dimf ,mn,nn,rf ,rdp,rdp2, 
:       nf ,mt,nf2,mt2 

PARAMETER   (DP  =  512) 

real  pi , t in(dp),ain(dp),rawa(dp),rwdf (dp*2,dp), rswdf (dp*2,dp), 
:     rawt(dp),dtsun,delt,dt  ,asuin,meanv,mt  ime.deM  , del2, 
:     const,  t  ,si_m,sunt>,val  ,sval  ,  udf  (dp*2,dp),coef  ,df  ,area 

complex  s(dp*2),dum,c(dp*5),dum3,dum2 

character*25  inname 

common/ in/  mvopt,redopt,dimt,dimf ,irm,nnf rf ,rdp,rdp2, 
:       nf ,mt,nf2,mt2 

common/ rl/  pi ,  t  i  n,  a  i  n, rawa , rwdf , rswdf , 
:     rawt,dtsum,delt,dt,asumfmeanv,mtime,del1,del2, 
:     const , t , sum, sumb, va I , sva I , wdf , coef , df , area 

common/comp/  s,dum,c,dum3,dum2 

dp2  =  dp*2 


df=1./(4.*dp*dt) 

rdp=dp/mm 
rdp2=dp2/nn 

nf2=nf*2 
mt2=mt*2 

val=1./((2.*pi)**2*nf*mt) 

do  20  j=-mt2,mt2 
do  10  i=-nf2,nf2 

coef=-((j*j)/(2.*mt*mt)+(i*i)/(2.*nf*nf)) 
hg( i , j )=val*exp(coef ) 
10     continue 
20   continue 

do  100  j=1,rdp 
do  100  i=1,rdp2 
100  rswdf (i,j)=0.0 


ii=0 

DO  500  N=1,DP2,nn 
jj  ■  0 
ii  =  ii+1 
DO  450  M  =  1,DP,mm 
jj  =  JJ+1 
DO  400  K=m-MT,m+MT 

DO  350  L  =  n-nf(n+nf 

IF  (L.LT.1)  THEN 

RSWDF(ii.jj)  =  RSWDF(ii,jj)  +  0 
ELSEIF  (K.LT.1)  THEN 
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RSUDFOi.jj)   =  RSUDFO'i.jj)   +  0 
ELSEIF   (L.GT.DP2)   THEN 

RSUDFO'i.jj)   =  RSUDF(ii,jj)  ♦   0 
ELSEIF   (K.GT.DP)   THEN 

RSUDFO'i.jj)  =  RSWDF(ii,jj)  +  0 
ELSE 

RSUDFO'i.jj)   =  RSUDF(ii,jj)+WDF(L,K)*HG(L-N,K-M) 
END  IF 


350 

CONTINUE 

400 

CONTINUE 

RSUDFO"  i ,  j  j  )=RSUDF(  i  i ,  j  j  ) 

450 

CONTINUE 

500 

CONTINUE 

END 
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SUBROUTINE  UIGH3 


*  This  subroutine  calculates  the  WDF  of  the  signal   * 

*  * 
********************************************************* 
* 

* 
* 

integer  dp,dp2,mvopt,redopt,dimt,dimf ,mm,m,rf , rdp,rdp2, 
:       nf fmt,nf2,mt2 
PARAMETER   (DP  =  512) 

real  pi , tin(dp),ain(dp),rawa(dp),rwdf (dp*2, dp), rswdf (dp*2,dp), 
:     rawt(dp),dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
:     const , t , sum, sumb, va I , sva I , wdf (dp*2 , dp ) , coef , df .area 
complex  s(dp*2),dum,c(dp*5),dum3,dum2 
character*25  inrvame 
* 

common/ in/  mvopt,redopt,dimt,dimf .mm.nn.rf ,rdp,rdp2, 
:       nf ,mt,nf2,mt2 

common/rl/  pi, tin,ain,rawa,rwdf, rswdf , 
:     rawt,dtsum,delt,dt,asum,meanv,mtime,del1,del2, 
:     const , t , sum, sumb, va I , sva I , wdf , coef , df , area 

common/ comp/  s,dum,c,dum3,dum2 

dp2  =  dp*2 


do  100  j  =  1  ,  dp 

*******auto  correlation  of  the  signal  corr3 
coef  =  2.0  *  dt 
do  300  i  =  1  ,  dp+1 

if(j.ge.i)  then 
dum  =  s(j-i*1) 
else 

dum  ■  cmplx(0.,0.) 
endif 

c(i)  =  coef  *  (s( j+i-1)*conjg(dum)) 

if (i .ne. j.and.i .ne.dp+1)  then 

c(dp2-i+2)  =  conjg(c(i)) 
endif 

300     cont  i  nue 

call  fft3 

do  200  i=1,dp2 

wdf(i,j)=real(c(i)) 
200       continue 

1 00   cont  i  nue 

END 
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C.        PLOTTING  PROGRAMS 

Listed  below  are  the  two  computer  programs  that  were  used  to  produce  the  Figures  in 
this  work. 

1.       2-D  Plotting  Routine 


PROGRAM  SIGNALPLT 
* 

********************************************************* 
* 

*  THIS  PROGRAM  USES  THE  GRAPHIC  PACKAGE  DISSPLA  TO 

*  PLOT  THE  A  2-D  FUNCTION. 
* 
********************************************************* 


***---DECLARING  VARIABLES---*** 

REAL  XARAY ( 1 000 ),YARAY( 1000) 

REAL  X,Y,F,XSTART,XEND,YSTART,YEND,XSTEP 

INTEGER  IPROPT,NPNTS,IX,IY 

*** INITIALIZING  VARIABLES  AND  SETTING  GRAPH  OPTIONS---*** 


— se 

WRITE(*,* 
WRITE(*,* 
READ(*,*) 
WRITE(*,* 
READ(*,*) 

—  se 
WRITE(*,* 
WRITE(*,* 
READ(*,*) 

—  se 
WRITE(*,* 
WRITE(*f* 
READ(*,*) 
WRITE(*,* 
READ(*,*) 

—  se 
WRITE(*f* 
WRITE(*,* 
READ(*,*) 

— se 
WRITE(*,* 
WRITE(*,* 
READ(*,*) 

— se 
WRITE(*,* 
WR1TE(*,* 
READ<*,*) 
WRITE(*,* 
READ<*,*> 

—  se 

WRITE(*.* 
WRITE(*,* 
WRITE(*,* 
READ(*,*), 


option  for 
•DECIDE  THE 
■INPUT  THE 

XSTART 
•INPUT  THE 

XEND 
desired  x  s 
•DECIDE  THE 
'INPUT  THE 
XSTEP 
option  for 
'DECIDE  THE 
'INPUT  THE 

YSTART 
'INPUT  THE 

YEND 
option  for 
'DECIDE  THE 
'INPUT  THE 
NPNTS 

option  for 
'DECIDE  THE 
'INPUT  THE 
I  MARK 

option  for 
'DECIDE  THE 
'INPUT  THE 

IX 
'INPUT  THE 

IY 


x  starting  &  ending  value--- 

START1NG  AND  ENDING  VALUES  OF  X1 
STARTING  VALUE  OF  X1 

ENDING  VALUE  OF  X' 

tep  size — 
X  STEP  SIZE' 
X  STEP  SIZE' 

y  starting  &  ending  value — 

STARTING  AND  ENDING  VALUES  OF  Y1 
STARTING  VALUE  OF  Y' 

ENDING  VALUE  OF  Y' 

number  of  points  to  plot  — 

NUMBER  OF  POINTS  TO  PLOT' 
NUMBER  OF  POINTS' 

line  style  (IMARK)--- 
IMARK  STYLE  DESIRED' 
IMARK  STYLE' 

grid  style — 

GRID  STYLE  DESIRED' 
X  GRID  STYLE' 

Y  GRID  STYLE' 


option  for  output  device — 

'YOU  WILL  NOW  DECIDE  WHERE  TO  SEND  THE  OUTPUT* 

'ENTER  0  FOR  OUTPUT  TO  THE  SCREEN' 

'ENTER  1  FOR  OUTPUT  TO  THE  LASER' 

IPROPT 
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* 

C  CALL  PDEV('LN03',IEER) 

IF  (IPROPT.EQ.O)  THEN 
CALL  PGPX 
ELSE 

CALL  LN03I 
END  IF 


"FUNCTION  EXECUTION  TO  DETERMINE  PLOTTING  POINTS** 

OPEN  (20,FILE='SIN400256A2.DAT',STATUS='OLD,) 
CALL  F(XARAY, YARAY, NPNTS, XSTART, XEND) 


***---LEVEL  ONE  COMMANDS---*** 
CALL  HWROT( 'COMIC') 
CALL  PAGE(11.0,8.5) 
CALL  NOBRDR 
CALL  AREA2D(9. 0,6.0) 

***-- -LEVEL  TWO  COMMANDS---*** 

CALL  SUISSB 

CALL  HEADINC400  HZ  SINE  WAVES' ,  100,1 .5,2) 

CALL  HEADIN( 'SIGNAL  LENGTH  =  0.128  SEC  /  AMP  =  2.0$' , 100,1 .5,2) 

CALL  SWISSM 

CALL  XNAMECTIME   (sec)$',100) 

CALL   YNAMEC AMPLITUDES ',100) 

CALL  GRAF(XSTART,XSTEP,XEND,YSTART, 'SCALE', YEND) 

***---LEVEL   THREE   COMMANDS---*** 

CALL  CURVE(XARAY,YARAY,NPNTS,IMARK) 

CALL  FRAME 

CALL  GRID(IX.IY) 

***-- -CLOSING  COMMANDS---*** 

CALL  ENDPL(O) 
CALL  DONEPL 
END 


*  ---FUNCTION  SUBROUTINE  TO  DETERMINE  PLOTTING  POINTS--- 

* 

SUBROUT I NE  F ( XARAY , YARAY , NPNTS , XSTART , XEND ) 

REAL  XARAY(NPNTS),YARAY(NPNTS), XSTART, XEND, STEP, XTEMP 
INTEGER  NPNTS 

*  ---READ  IN  THE  XARAY  AND  YARAY  VALUES--- 
DO  100  I  =  1, NPNTS 

READ(20,*)  XARAY(I),  YARAY(I) 
100    CONTINUE 
END 
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2.       3-D  Plotting  Routine 


PROGRAM  UDFPLTVER3 

* 
********************************************************* 

* 

*  THIS  PROGRAM  USES  THE  GRAPHIC  PACKAGE  DISSPLA  TO 

*  PLOT  A  3-D  FUNCTION. 

* 
********************************************************* 

* 

*  "'---DECLARING  VARIABLES---*** 

* 

REAL  RWDF(32768) 

INTEGER  IPROPT,IXPTS,IYPTS 

* 

*  *** — INITIALIZING  VARIABLES  AND  SETTING  GRAPH  OPTIONS---*** 

OPEN(15,FILE='RSWDF.OUT',STATUS='OLD') 

* 

DO  100  I  =  1,32768 
READ(15,*)  RUDF(I) 
100      CONTINUE 


call  pdev('ln03',  ieer) 

♦♦FUNCTION  EXECUTION  TO  DETERMINE  PLOT** 

***---LEVEL  ONE  COMMANDS---*** 

CALL  PAGE(8.5,11.0) 
CALL  AREA2D(8. 0,7.0) 

***--- LEVEL  TWO  COMMANDS---*** 

CALL  VOLM3D(6.,6.,4.) 

CALL  HEADINC200  &  600  Hz  SINE  WAVES' ,  100, -1 .5,2) 

CALL  HEAD  INC  FREQUENCY  RESOLUTION  ;  5  X  35$' , 100, -1 .5,2) 

CALL  X3NAMEC FREQUENCY* ',100) 

CALL  Y3NAMEC  ',1) 

CALL  Z3NAMEC AMPLITUDES ',100) 

CALL  VUANGL(-90.,0.,30.) 

CALL  GRAF3D(0.,1250./5., 1250.0, 
:  0.,. 2048/2., 0.2048, 

:  0., 'SCALE', .01) 

***-- -LEVEL  THREE  COMMANDS---*** 

CALL  SURMAT(RWDF, 1,256, 1,128,0) 
CALL  FRAME 

***... CLOSING  COMMANDS---*** 

CALL  ENDPL(O) 
CALL  DONEPL 

END 
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